!#############################################################
 module share
! Shared variables, constants and a few generic subroutines
!#############################################################
   implicit none
   integer :: rgx,rgz,nx,nz,Nmk,timin,timax,stp_ascii,write_mks,write_mm_matrix,write_e_matrix,                    &
              debug,solver,mord,minfo,ti,nli,maxit,nlv,conv_test,cpr,tala,ala,therm,markers,mkinterp,              &
	      jph_400,jph_670,restart,nphases,edepth,visc_law,overwr_visc,overwr_visc_mod,                         &
	      alpha_law,kappa_law,alpha_temp,kappa_temp,                                                           &
              icd,itd,ihd,ivd,decayQ,core_cooling,uTconst,uBconst,vLconst,vRconst,bcB_T,melt,phase_dep_ap,         &
              bcB_flow,bcT_flow,bcL_flow,bcR_flow

   integer, dimension(:), allocatable :: bcB,bcT,bcL,bcR

   integer, dimension(:,:), allocatable :: visc_mech

   real (kind=8) :: dx,dz,dimx,dimz,Cnr,time_ascii,relax_fac,ext_fac,                                              &
                    dt,init_dt,time_max,nlitol,Ra,RaC,Di,RaHval,Hdimz,vrmsold,dph,                                 & 
                    drho_400,drho_670,drho_ppv,Tph_400,Tph_670,T_ppv,kpvppv,eolsp,esppv,epvppv,Ft,Fz,              &  
                    ae,av,ae_ds,av_ds,stress_exp,epsref,prefac,prefac_ds,Tref,Zref,sigmay,sigmayz,yield_depth,etamax,etamin, & 
                    decay_lambda,cooling_factor,age,Tbott,Ttop,Tsurf,qbott,Tsurf_d,Tcmb_d,deltaS,lidfac,mixdepth,  &
                    alpha0,alpha1,alpha2,alpha3,kap0,kap1,kap2

   real (kind=8), dimension(:,:), allocatable :: u,v,p,dudx,dvdz,dudz,dvdx,Phi,alpha,kap,urms,RaH,RaH_old,         &
                                                 uold,vold,uextr,vextr,Told,Tnew,Tpr,Gam,Gam_old,eta,eta_old,      &
                                                 drho,depl,depl_old,Tsol,Tsol_old
 
   real (kind=8), dimension(:), allocatable :: uT,uB,vL,vR,Tsol0,Tliq,                                 &
                                               phase_Rb,phase_z,phase_Cl,phase_T,phase_eta,            &
                                               phase_alpha0,phase_alpha1,phase_alpha2,phase_alpha3,    &
                                               phase_kap0,phase_kap1,phase_kap2,                       &
                                               phase_ae,phase_av,phase_ae_ds,phase_av_ds,phase_prefac, &
                                               phase_prefac_ds

   real (kind=8), dimension(:,:,:), allocatable :: phase_f

   integer :: comp_mk, heats_mk, visc_mk, depl_mk

   real (kind=8), dimension(:), allocatable :: xmk,zmk,xgrid,zgrid,delx,delz,mktype,Gmk,dGmk,RaHmk,etamk,deplmk,umk,vmk,Tmk

   real (kind=8), parameter :: pi=3.14159265358979d0

   character(len=100) :: outdir

!*******************
  contains
!*******************

!================================================================
 subroutine read_input_file
! Read input variables from yacc.inp
!================================================================
   integer :: i
   character (len=100) :: readval
   character (len=1)  :: ch
!=========================================

!--------------------------------------------
! Output stuff
!--------------------------------------------
   call readinput("output_dir",readval)
   outdir=readval

   call readinput("debug",readval)
   read(readval,*) debug

   call readinput("stp_ascii",readval)
   read(readval,*) stp_ascii

   call readinput("time_ascii",readval)
   read(readval,*) time_ascii

   call readinput("write_mks",readval)
   read(readval,*) write_mks

!---------------------------------------------
! Boundary conditions
!---------------------------------------------
! Energy (bottom)
   call readinput("bottom_bc_temp",readval)
   read(readval,*) bcB_T

! Momentum top
   call readinput("top_bc_flow",readval)
   read(readval,*) bcT_flow

   if(bcT_flow.eq.1) then 
     call readinput("top_u_const",readval)
     read(readval,*) uTconst
   endif

! Momentum bottom
   call readinput("bottom_bc_flow",readval)
   read(readval,*) bcB_flow

   if(bcB_flow.eq.1) then 
     call readinput("bottom_u_const",readval)
     read(readval,*) uBconst
   endif

! Momentum left
   call readinput("left_bc_flow",readval)
   read(readval,*) bcL_flow

   if(bcL_flow.eq.1) then 
     call readinput("left_v_const",readval)
     read(readval,*) vLconst
   endif

! Momentum right
   call readinput("right_bc_flow",readval)
   read(readval,*) bcR_flow

   if(bcR_flow.eq.1) then 
     call readinput("right_v_const",readval)
     read(readval,*) vRconst
   endif
   
!--------------------
! Time evolution
!--------------------
   call readinput("restart",readval)
   read(readval,*) restart

   call readinput("nl_viscosity",readval)
   read(readval,*) nlv

   if(nlv.eq.1) then
     call readinput("conv_test",readval)
     read(readval,*) conv_test

     call readinput("relax_factor",readval)
     read(readval,*) relax_fac

     call readinput("max_iterations",readval)
     read(readval,*) maxit
   endif

   call readinput("nlitol",readval)
   read(readval,*) nlitol

   call readinput("energy",readval)
   read(readval,*) therm

   call readinput("time_max",readval)
   read(readval,*) time_max

   call readinput("timax",readval)
   read(readval,*) timax

   call readinput("courant",readval)
   read(readval,*) Cnr

   call readinput("time_extrap_factor",readval)
   read(readval,*) ext_fac

   call readinput("time_step",readval)
   read(readval,*) dt

   call readinput("init_time_step",readval)
   read(readval,*) init_dt

   call readinput("decayQ",readval)
   read(readval,*) decayQ

   call readinput("decay_lambda",readval)
   read(readval,*) decay_lambda

   call readinput("core_cooling",readval)
   read(readval,*) core_cooling

   if(core_cooling.eq.1) then
     call readinput("cooling_factor",readval)
     read(readval,*) cooling_factor
   endif

!--------------------------
! Fields on markers
!--------------------------
   if(markers.eq.1) then
     call readinput("comp_mk",readval)
     read(readval,*) comp_mk   

     call readinput("heats_mk",readval)
     read(readval,*) heats_mk 

     call readinput("visc_mk",readval)
     read(readval,*) visc_mk 
   endif

!---------------------------------------------
! Basic equations and non-dimensional groups
!---------------------------------------------
   call readinput("TALA",readval)
   read(readval,*) tala      ! TALA equations

   call readinput("ALA",readval)
   read(readval,*) ala       ! ALA equations

   if(tala.eq.1.or.ala.eq.1) then
     cpr=1
   endif
   
   call readinput("Ra",readval)
   read(readval,*) Ra            ! Rayleigh number

   call readinput("RaC",readval)
   read(readval,*) RaC           ! Compositional Rayleigh number

   call readinput("RaH",readval)
   read(readval,*) RaHval        ! Rayleigh number for internal heating

   call readinput("Di",readval)
   read(readval,*) Di            ! dissipation number

!------------------------
! Viscosity parameters
!------------------------
   call readinput("viscosity_law",readval)
   read(readval,*) visc_law     ! Viscosity law (1: FK, 2: Arr)

   call readinput("overwrite_viscosity",readval)
   read(readval,*) overwr_visc  ! (Partially) prescribed viscosity

   if(overwr_visc.eq.1) then
     call readinput("overwrite_viscosity_model",readval)
     read(readval,*) overwr_visc_mod      
   endif

   call readinput("T_surface",readval)
   read(readval,*) Tsurf        ! non-dimensional surface temperature 

   call readinput("max_viscosity",readval)
   read(readval,*) etamax  ! maximum viscosity

   call readinput("min_viscosity",readval)
   read(readval,*) etamin  ! minimum viscosity

   call readinput("ae",readval)
   read(readval,*) ae           ! Activation energy for diffusion creep

   call readinput("av",readval) 
   read(readval,*) av           ! Activation volume for diffusion creep

   call readinput("ae_ds",readval)
   read(readval,*) ae_ds         ! Activation energy for dislocation creep

   call readinput("av_ds",readval) 
   read(readval,*) av_ds         ! Activation volume for dislocation creep

   call readinput("stress_exp",readval) 
   read(readval,*) stress_exp    ! Stress exponent for dislocation creep

   call readinput("prefac",readval) 
   read(readval,*) prefac        ! prefactor for diffusion creep

   call readinput("prefac_ds",readval) 
   read(readval,*) prefac_ds     ! prefactor for dislocation creep 

   call readinput("eT",readval)
   read(readval,*) Ft           ! Viscosity contrast due to T fro FK law

   call readinput("ez",readval)
   read(readval,*) Fz           ! Viscosity contrast due to p fro FK law 

   call readinput("sigma_yield",readval)
   read(readval,*) sigmay       ! yield stress

   call readinput("sigma_yield_z",readval)
   read(readval,*) sigmayz       ! yield stress gradient

   call readinput("yield_depth",readval)
   read(readval,*) yield_depth   ! max depth where yielding is considered

   call readinput("Tref",readval)
   read(readval,*) Tref         ! Reference temperature

   call readinput("Zref",readval)
   read(readval,*) Zref         ! Reference depth

   call readinput("epsref",readval)
   read(readval,*) epsref       ! Reference strain rate

   call readinput("staglid_factor",readval)
   read(readval,*) lidfac       ! Factor to calculate stagnant lid

   call readinput("mixdepth",readval)
   read(readval,*) mixdepth     ! Depth below which the mixed fraction is calculated

!------------------------------------------
! Melting variables
!------------------------------------------
   call readinput("melt",readval)
   read(readval,*) melt

   if(melt.eq.1) then
     call readinput("depl_mk",readval)
     read(readval,*) depl_mk 

     call readinput("deltaS",readval)
     read(readval,*) deltaS
   endif

!------------------------------------------
! Phase transitions
!------------------------------------------
   call readinput("phases",readval)
   read(readval,*) nphases       ! number of phases

   if(nphases.ne.0) then

     call readinput("phase_loop",readval)
     read(readval,*) dph       ! amplitude of the phase transition boundary  

     allocate(phase_Rb(nphases),phase_z(nphases),phase_Cl(nphases),phase_T(nphases),phase_eta(nphases), &
              phase_f(nphases,0:nx+1,0:nz+1))

     do i=1,nphases
   
       write (ch, '(I1)') i

       call readinput("phase_Rb_"//ch,readval)
       read(readval,*) phase_Rb(i)

       call readinput("phase_z_"//ch,readval)
       read(readval,*) phase_z(i)

       call readinput("phase_Cl_"//ch,readval)
       read(readval,*) phase_Cl(i)

       call readinput("phase_T_"//ch,readval)
       read(readval,*) phase_T(i)

       call readinput("phase_eta_"//ch,readval)
       read(readval,*) phase_eta(i)

     end do

!   Phase dependent activation parameters  
     call readinput("phase_dep_ap",readval)
     read(readval,*) phase_dep_ap       

     if(phase_dep_ap.eq.1) then
       allocate(phase_ae(nphases+1),phase_av(nphases+1),phase_ae_ds(nphases+1),phase_av_ds(nphases+1), &
                phase_prefac(nphases+1),phase_prefac_ds(nphases+1))

       do i=1,nphases+1
   
         write (ch, '(I1)') i

         call readinput("phase_ae_"//ch,readval)
         read(readval,*) phase_ae(i)

         call readinput("phase_av_"//ch,readval)
         read(readval,*) phase_av(i)

         call readinput("phase_ae_ds_"//ch,readval)
         read(readval,*) phase_ae_ds(i)

         call readinput("phase_av_ds_"//ch,readval)
         read(readval,*) phase_av_ds(i)

         call readinput("phase_prefac_"//ch,readval)
         read(readval,*) phase_prefac(i)

         call readinput("phase_prefac_ds_"//ch,readval)
         read(readval,*) phase_prefac_ds(i)

       end do

     endif

   endif

!------------------------------------------
! Thermal expansivity
!------------------------------------------
   call readinput("alpha_law",readval)
   read(readval,*) alpha_law   

   call readinput("alpha_temp",readval)
   read(readval,*) alpha_temp  

! pT-dependent expansivity (no phase change)
   if(alpha_law.eq.1) then
     call readinput("alpha0",readval)
     read(readval,*) alpha0  
     call readinput("alpha1",readval)
     read(readval,*) alpha1  
     call readinput("alpha2",readval)
     read(readval,*) alpha2  
     call readinput("alpha3",readval)
     read(readval,*) alpha3
   endif

! pT-dependent expansivity (different coefficients for each phase)
   if(nphases.ne.0.and.alpha_law.eq.2) then

     allocate(phase_alpha0(nphases+1),phase_alpha1(nphases+1),phase_alpha2(nphases+1),phase_alpha3(nphases+1))

     do i=1,nphases+1
   
       write (ch, '(I1)') i

       call readinput("phase_alpha0_"//ch,readval)
       read(readval,*) phase_alpha0(i)

       call readinput("phase_alpha1_"//ch,readval)
       read(readval,*) phase_alpha1(i)

       call readinput("phase_alpha2_"//ch,readval)
       read(readval,*) phase_alpha2(i)

       call readinput("phase_alpha3_"//ch,readval)
       read(readval,*) phase_alpha3(i)

     end do

   endif

!------------------------------------------
! Thermal conductivity
!------------------------------------------
   call readinput("k_law",readval)
   read(readval,*) kappa_law  
   
   call readinput("k_temp",readval)
   read(readval,*) kappa_temp  

! pT-dependent conductivity (no phase change)
   if(kappa_law.eq.1) then

     call readinput("k0",readval)
     read(readval,*) kap0  
     call readinput("k1",readval)
     read(readval,*) kap1  
     call readinput("k2",readval)
     read(readval,*) kap2  

! pT-dependent conductivity (different coefficients for each phase)
   elseif(nphases.ne.0.and.kappa_law.eq.2) then

     allocate(phase_kap0(nphases+1),phase_kap1(nphases+1),phase_kap2(nphases+1))

     do i=1,nphases+1
   
       write (ch, '(I1)') i

       call readinput("phase_k0_"//ch,readval)
       read(readval,*) phase_kap0(i)

       call readinput("phase_k1_"//ch,readval)
       read(readval,*) phase_kap1(i)

       call readinput("phase_k2_"//ch,readval)
       read(readval,*) phase_kap2(i)

     end do
     
   endif

!------------
! Solver
!------------
   call readinput("solver",readval)
   read(readval,*) solver

   call readinput("matrix_ord",readval)
   read(readval,*) mord

   call readinput("matrix_info",readval)
   read(readval,*) minfo

   call readinput("write_mm_matrix",readval)
   read(readval,*) write_mm_matrix

   call readinput("write_e_matrix",readval)
   read(readval,*) write_e_matrix


 end subroutine read_input_file

 
!=============================================================================
 subroutine readinput(string,value)
! Search the input file and assigns a value to a selected string variable
!=============================================================================
   implicit none
   integer :: k,j
   character(len=*) :: string,value
   character(len=100) :: s,v
!=============================================================================

!---------------------------
! Read from input file
!---------------------------
   k=0
   open(1,file="yacc.inp")
   
   do while(.true.)

     read(1,*,end=99) s
     if(s(1:1).ne.'#') then
       backspace(1)
       read(1,*,end=99) s,v
     endif
     
     if (s.eq.string) then
       value=v
       k=1
     endif
     
   end do
99 continue
   close(1)

!----------------------------------------------------------------
! If a parameter is not found, it is read from the default file
!----------------------------------------------------------------
   j=0
   if(k.ne.1) then

     open(2,file="yacc.def")
     do while(.true.)

       read(2,*,end=88) s,v
     
       if (s.eq.string) then
         value=v
         j=1
         print *, 'Variable ',string,' missing in yacc.inp'
         write(*,'(A13,3X,A8,A22)') 'Default value',v, 'assigned from yacc.def'
       endif
     
     end do
88   continue
     close(2)
   endif

!----------------------------------------------------------------
! If also yacc.def misses the parameter, stop
!----------------------------------------------------------------
   if(k.ne.1.and.j.ne.1) then
     print *, 'Variable ',string,' missing in yacc.inp and yacc.def => STOP'
     stop
   endif

 
 end subroutine readinput


!===========================================
 function ran(idum)
! Park-Miller random number generator
! see Numerical Recipes Ch. 7
!===========================================
   implicit none
   integer, parameter :: k4b=selected_int_kind(9)
   integer(k4b), intent(inout) :: idum
   real :: ran
   integer(k4b), parameter :: ia=16807,im=2147483647,iq=127773,ir=2836
   real, save :: am
   integer(k4b), save :: ix=-1,iy=-1,k

   if (idum <= 0 .or. iy < 0) then
        am=nearest(1.0,-1.0)/im
        iy=ior(ieor(888889999,abs(idum)),1)
        ix=ieor(777755555,abs(idum))
        idum=abs(idum)+1
   end if
   ix=ieor(ix,ishft(ix,13))
   ix=ieor(ix,ishft(ix,-17))
   ix=ieor(ix,ishft(ix,5))
   k=iy/iq
   iy=ia*(iy-k*iq)-ir*k
   if (iy < 0) iy=iy+im
   ran=am*ior(iand(im,ieor(ix,iy)),1)

 end function ran

!========================================================================
 function bisec(xx,x)
! Bisection algorithm to locate the position of 'x' on an 
! ascending grid contained in 'xx'
! On output the index 'jl' is returned such that xx(jl) <= x <= xx(jl+1)
! 
! Input:
! xx => real, array(0:). Either x- or z-grid
! x => real. x- or z- coord of the point to be located
! Output: 
! bisec => integer. Left index of the bounding cell
!========================================================================
   implicit none
   integer :: bisec,jl,ju,jm,n
   real(kind=8), dimension(0:), intent(in) :: xx
   real(kind=8), intent(in) :: x

   n=size(xx)-2
   jl=0
   ju=n+1
   do
     if(ju-jl.le.1) exit
     jm=(ju+jl)/2
     if(x.ge.xx(jm)) then
       jl=jm
     else
       ju=jm
     end if
   end do
   
   if(x.eq.xx(0)) then
     bisec=0
   elseif(x.eq.xx(n+1)) then
     bisec=n+1
   else
     bisec=jl
   end if
   
 end function bisec

!=============================================
 subroutine errorf(x,err)
! Compute error function erf(x)
! Input: 
! x => real. Argument of the erf(x)
! Output:  
! err => real. Error function of x
!=============================================
   implicit none
   integer :: k
   real (kind=8) :: eps,pi,x,x2,err,er,r,c0
!=============================================
     
   eps=1.0d-15
   pi=3.141592653589793d0
   x2=x*x
   if (dabs(x).lt.3.5d0) then
      er=1.0d0
      r=1.0d0
      do k=1,50
         r=r*x2/(k+0.5d0)
         er=er+r
         if (dabs(r).le.dabs(er)*eps) exit
      end do
      c0=2.0d0/dsqrt(pi)*x*dexp(-x2)
      err=c0*er
   else
      er=1.0d0
      r=1.0d0
      do k=1,12
         r=-r*(k-0.5d0)/x2
         er=er+r
      end do  
      c0=dexp(-x2)/(dabs(x)*dsqrt(pi))
      err=1.0d0-c0*er
      if (x.lt.0.0) err=-err
   endif
   return

 end subroutine errorf


 
 end module share
!###############################################

